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I.  INTRODUCTION 


The  purpose  of  this  effort  was  to  modify  an  existing  computer  code 
(Ref.  1)  for  more  efficient  computation  of  viscous  flowfields  in  and  about 
hollow  projectile  shapes,  (See  Fig.  1.)  The  existing  code  solves  the  axi- 
symmetric  thin-layer  equations  on  domains  consisting  of  a  single  element 
or  module.  A  typical  grid  for  a  problem  of  this  type  is  shown  in  Fig.  2. 

Grids  of  this  type  suffer  from  several  problems.  First  of  all,  skewness 
problems  in  such  a  grid  prevent  accurate  solutions  and  can  even  lead  to 
instabilities.  Secondly,  computed  inefficiency  exists  due  to  the  concentration 
of  points  where  they  are  not  needed.  Finally,  this  type  of  grid  results  in 
severe  and  undesirable  gradients  in  geometric  parameters  near  the  inflow 
boundary  intersection  with  the  projectile  axis.  Of  course  these  grid  prob¬ 
lems  can  be  quite  detrimental  to  the  quality  of  the  computed  solution  and 
can  result  in  excessive  CPU  time  requirements. 

The  method  chosen  to  eliminate  these  grid  problems  was  to  develop  a 
modular  grid  system.  The  schematic  of  such  a  system  is  shown  in  Fig.  3.  A 
grid  system  of  this  type  has  several  attractive  features.  First,  points  are 
concentrated  where  they  are  needed.  For  example,  sufficient  refinement  is 
obtained  in  the  boundary -layer  regions  without  forced  refinement  where  it  is 
not  needed.  Also, the  number  of  vertical  points  in  Region  2  is  independent 
of  that  in  Region  3.  This  prevents  the  wasted  computation  on  the  inside  of 
the  projectile  that  existed  for  the  original  grid.  Secondly,  the  geometric 
parameters  are  well  behaved  in  the  corner  where  the  inflow  boundary  inter¬ 
sects  the  projectile  centerline.  In  addition,  this  type  of  grid  system 
virtually  eliminates  skewness  problems. 
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There  are  some  disadvantages  of  this  type  of  grid  system,  also.  First, 
it  requires  grid  generation,  flow  solution,  output,  etc.,  for  each  distinct 
region  rather  than  just  one  domain.  In  addition,  the  interfaces  between 
the  distinct  regions  of  the  grid  require  special  consideration  since  the 
separate  grids  need  not  have  continuous  mappings  across  these  interfaces. 

For  example,  between  Regions  1  and  2  the  grid  suddenly  changes  from  a  very 
fine  mesh  to  a  coarse  one.  As  a  result  of  these  special  considerations,  the 
coding  logic  for  performing  the  block-tridiagonal  inversions  through  these 
interfaces  becomes  very  complex  and  is  compounded  further  by  the  fact  that 
the  strong-conservation- law  form  of  the  governing  fluid  equations  is 
employed.  This  fact  will  become  obvious  in  the  development  of  the  required 
interface  operators. 

The  fundamental  notion  employed  in  this  work  is  that  second-order 
accurate  central-difference  approximations  for  first  and  second  derivatives 
on  a  uniformly  spaced  grid  in  computational  space  may  be  replaced  by  second- 
order  accurate  approximations  on  an  unequally  spaced  grid  in  computational 
space  without  compromising  the  tridiagonal  structure  of  the  difference 
operator.  This  is  demonstrated  in  the  following  section. 


II.  CONCEPT  TESTING 


The  notion  of  employing  finite-difference  approximations  to  an  unequally 
spaced  grid  in  computational  space  is  tested  on  a  one -dimensional  model  prob¬ 


lem  governed  by  the  nonlinear  viscous  Burgers '  equation 
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subject  to  the  boundary  conditions 
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u(0,t)  =  1 


u(l,t)  =  0 


and  the  initial  conditions 


1  ,  x  =  0 


u(x,0) 


0  ,  0  <  x  <  1 


The  exact  steady-state  solution  to  this  problem  is 


u  =  u  tanh|-jjj  (1  -  x)  j 


u  tanh 


where  u  solves  the  transcendental  equation 


A  -  A 

u  -  1  -u/y 

A  «  “  G 

U  +  1 


This  problem  was  chosen  for  several  reasons.  First,  the  exact  solution  is 
known;  thus,  a  precise  measure  of  the  error  generated  by  the  finite-difference 
scheme  is  known.  Second,  for  sufficiently  small  values  of  JJ,  a  steep  gradi¬ 
ent  in  the  solution  exists  near  the  right  boundary  of  the  domain.  (See  Fig.  4.) 
This  behavior  is  similar  to  that  which  exists  in  a  boundary-layer  region. 

Thus  the  domain  can  be  partitioned  into  two  grids:  a  fine  mesh  near  the 
right  boundary  and  a  coarse  mesh  elsewhere.  Since  y  can  be  adjusted  to 
position  the  steep  gradient  either  completely  within  the  fine  mesh  or  so  that 
it  crosses  the  intermesh  boundary,  the  required  intermesh  operators  can  be 
severely  tested  with  this  problem. 

Equation  (1)  must  be  transformed  to  computational  space  according  to 
the  mapping  T  =  t,t  £  =  £(x).  Note  that  actually  two  mappings  are  employed, 
one  for  the  fine  grid  and  one  for  the  coarse  grid.  The  mapped  equation  is 
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Consistent  with  the  tubular  projectile  code,  the  Beam-Warming  implicit 


numerical  integration  scheme  in  delta  form  (Ref.  2)  is  used  to  solve  this 
equation.  The  scheme  is  expressed  as 
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where  3  =  1  yields  Euler  implicit  time  differencing  and  3=1/2  yields 
trapezoidal  time  differencing.  Conventional  central-difference  approxi¬ 
mations  are  used  interior  to  both  fine  and  coarse  meshes.  However ,  at  the 
intermesh  point  where  x  =  1  -  6  (see  Fig.  5) ,  a  special  formula  must  be 
applied  for  the  derivative  approximations.  This  is  derived  as  follows.  Con¬ 
sider  a  transformation  from  computational  space  back  to  physical  space  into 
an  arc  length  variable  such  that 

3_  _  3_6_  3 _ 

3  £  3?  34 

where  4  is  arc  length  along  the  £  direction.  Note  that  this  mapping  can  be 
performed  locally  at  x  =  1  -  6  from  either  the  fine  mesh  computational 
system  or  the  coarse  mesh  system.  The  3/94  derivative  is  then  approximated 
using  a  second-order  accurate  formula  for  an  unequally  spaced  grid.  The 
factor  3,4/9  £  simply  makes  the  difference  approximation  consistent  with  which¬ 
ever  grid  the  3/9  E,  is  meant  to  apply  to.  Of  particular  importance  is  the 
fact  that  this  procedure  will  work  on  a  multidimensional  problem  provided 
continuity  of  slope  across  interfaces  exists  for  all  grid  lines.  Second 
derivatives  in  computational  space  transform  back  to  arc  length  dependence 
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according  to 


3i_=  (Z±f  *L+  2k !_ 

3S'  a^2  a^2  34 

The  arc  length  derivatives  are  differenced  according  to 


32 ( •)  2  [ A ( ♦ )  V(»)  1 

a42  -  A  +  V  l  A  V  J 


3(')  =  1 

3-6  A  +  V 


I  ac)  *f  vo] 


(3) 

(4) 


where  A,V  refer  to  conventional  forward  and  backward  difference  operators, 
respectively,  and  A,V  alone  indicate  operation  on  the  arc  length  parameter, 

-6.  For  example,  A  =  A4  =  ~  4 .  The  factors  34/3  £  and  3  4/3  £  are 

approximated  by  one-sided  finite  differences  using  the  arc  length  data  on 
whichever  side  is  appropriate. 

The  result  of  the  procedure  just  outlined,  insofar  as  the  tridiagonal 
algorithm  is  concerned,  is  simply  a  modification  to  the  algebraic  equation 
corresponding  to  the  intermesh  point.  This  requires  a  change  in  the  row 
of  the  tridiagonal  matrix  corresponding  to  this  point.  Also,  the  right-hand 
side  derivatives  at  this  point  must  be  differenced  according  to  the  procedures 
just  discussed.  Let  point  p  be  the  intermesh  point.  Then  p  +  1  is  in  the 
fine  mesh  and  p  -  1  is  in  the  coarse  mesh.  The  difference  scheme  represented 
by  Eq.  (2)  is  represented  at  this  point  as 

a  Au  „  +  b  Au  +  c  Au  „  =  K 
P  P-1  P  P  P  P+1  P 

where 
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where  superscript  *  indicates  that  the  quantity  is  evaluated  with  a  one¬ 
sided  difference  approximation  or  a  selected  side  of  the  intermesh  boundary. 
The  star  quantities  must  all  be  evaluated  on  the  same  side  for  consistency 
of  the  equation  although  which  side  is  selected  is  optional. 

This  intermesh  differencing  was  employed  and  various  test  cases  were 
run.  The  results  demonstrated  that  this  concept  is  a  valid  one.  Figure  6 
illustrates  the  steady-state  error  distribution  for  the  case  u  =  0.2,  6  =0.2 
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using  the  Euler  implicit  scheme  (6  =  1)  at  a  Courant  number  of  10.  Various 
values  of  y  and  6  were  used  and  the  results  didn't  differ  appreciably  from 
those  shown  so  long  as  y  0.1  and  no  adverse  clustering  was  used  in  either 
of  the  fine  or  coarse  mesh  regions.  Adverse  clustering  is  defined  as  follows. 
Define  A^  =  d  ~  <S)/Nc  and  A^  =  6/N^  where  Nc>  represent  the  number  of 
mesh  intervals  in  the  coarse  and  fine  mesh,  respectively.  Next,  form  the 
ratio  R  =  Af/Ac-  Clustering  in  either  region  which  causes 


R 


is  defined  as  adverse  clustering.  This  type  of  clustering  caused  the  error 
distribution  to  grow  but  was  not  of  concern  since  such  adverse  clustering  is 
not  used  in  practice.  Finally,  solutions  could  not  be  obtained  for  y  <  0.1 
no  matter  what  Courant  number  was  used.  This  was  true  even  for  the  limiting 
case  of  6=  0.5  and  Nc  =  Nv  which  corresponds  to  an  equally  spaced  mesh  over 
the  entire  domain.  Thus,  it  was  shown  that  this  problem  was  unrelated  to  the 
new  operators  for  unequally  spaced  grids.  It  is  suspected  that  this  problem 
is  a  result  of  using  central  differences  for  convective  terms  when  the 
effective  Reynolds  number  is  sufficiently  large  (y  sufficiently  small) . 

Summarizing  this  section,  the  concept  of  unequally  spaced  operators 
was  introduced  and  tested  on  a  model  equation.  The  results  of  this  test 
prove  the  feasibility  of  the  concept.  The  remaining  sections  expand  the  idea 
to  the  tubular  projectile  problem. 


13 


III.  GRID  GENERATION 


Figure  3  illustrates  the  new  modularized  grid  system  to  be  employed. 

The  actual  generation  of  the  grid  is  accomplished  with  a  hyperbolic  grid 
generation  scheme  as  described  in  Ref.  3.  The  first  step  is  to  generate  the 
fine  grid  for  Region  1  with  this  generator.  The  grid  for  Region  1  is 
marched  out  a  distance  sufficient  that  this  grid  contains  the  entire  boundary 
layer.  The  outer  boundary  of  Region  1  then  serves  as  a  portion  of  the  upper 
boundary  of  Region  3  and  the  lower  boundary  of  Region  2.  The  grid  line 
emanating  from  the  projectile  leading  edge  is  constrained  to  be  horizontal 
and  is  extended  to  the  inflow  boundary.  This  line  separates  Regions  2  and  3 
forward  of  the  fine  mesh  boundary.  The  grid  for  Region  2  is  then  generated 
from  the  starting  line  consisting  of  the  straight  section  forward  of  the  fine 
mesh  continued  by  the  upper  portion  of  the  outer  boundary  of  the  fine  mesh. 
This  starting  line  is  marched  up  to  a  distance  above  the  projectile  of  a 
few  chord  lengths.  The  grid  for  Region  3  is  obtained  in  a  similar  fashion 
except  that  the  starting  line  is  marched  down  until  every  point  corresponding 
to  the  last  £  =  constant  line  is  below  the  projectile  centerline.  Then  the 
points  along  each  £  =  constant  line  are  adjusted  along  this  curve  so  that 
the  last  £  =  constant  line  coincides  with  the  projectile  centerline.  The 
orthogonal  nature  of  this  grid  generator  provides  the  needed  slope  continuity 
condition  across  the  intermesh  boundaries.  It  should  be  noted  that  the  grid 
in  Region  3  is  not  quite  orthogonal  after  the  readjustment  process  to  fit  the 
to  the  projectile  centerline;  however,  it  is  nearly  so  in  practice.  An 
actual  grid  generated  with  this  scheme  is  shown  in  Fig.  7. 

Application  of  the  implicit  approximate-factored  integration  scheme  to 
a  grid  system  such  as  the  one  shown  in  Fig.  7  (schematic  in  Fig.  3)  is  far 


from  arbitrary.  There  exists  a  special  sequence  (not  unique)  of  operations 
which  taken  together  result  in  the  final  solution  on  the  existing  grid 
system.  The  structure  of  the  sequence  chosen  for  the  present  study  and 
some  motivation  for  this  choice  is  the  subject  of  the  next  section. 


IV,  INTEGRATION  SCHEME  TOPOLOGY 


The  integration  scheme  used  may  be  represented  as 

ft.JLAC  =  RHS 
12 

where  are  the  factored  finite-difference  operators  in  the  two  compu¬ 

tational  coordinate  directions.  The  order  in  which  the  directions  appear  is 
arbitrary  but  this  choice  can  be  made  to  result  in  easier  application  of 
boundary  conditions.  In  this  study  it  was  necessary  to  solve* =  RHS 
in  some  regions  and  =  RHS  in  others  for  simplicity.  In  general ,  the 

solution  algorithm  consists  of  four  major  steps: 

1.  Given  Q1 ,  compute  RHS. 

2.  Solve  £2^U*  =  RHS  for  U*  by  inverting  block  tridiagonal  system. 

3.  Solve  £22Aq  =  U*  for  Aq  by  inverting  block  tridiagonal  system. 

4.  Compute  new  Q  by  =  Qn  +  Aq. 

Steps  2  and  3  involve  the  £  and  £  direction  operators  which ,  as  mentioned , 

may  appear  in  either  order.  The  application  of  the  integration  scheme  to 

the  tubular  projectile  problem  with  the  present  modularized  grid  involves 

proper  sequencing  of  the  £  and  £  inversions.  This  is  best  described  by 

itemizing  each  step.  (Refer  to  Fig.  8.) 

1.  Perform  the  ^-inversion  on  interior  of  Region  3. 

Solve  =  RHS  for  U*  =  £2  Ao. 

£  C 
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2.  Perform  the  ^-inversion  on  the  interior  of  Region  2. 

Solve  fi^U*  «  RHS  for  U*  =  Q^Aq. 

3.  Perform  the  ^-inversion  in  Region  L  (Include  outer  boundary. ) 

Solve  fi>. U*  =  RHS  for  U*  =  fi^AQ. 

4.  Perform  the  ^-inversion  along  the  forward  cut  (leading  edge  grid 
line)  from  inflow  boundary  to  the  projectile  leading  edge  (line 
AC). 

5.  Perform  the  ^-inversion  along  the  trailing  edge  cut,  (See  line  DE 
of  Fig.  8.) 

6.  Perform  the  C-inversion  for  leading  edge  region  !a! 

(from  H  to  I).  (Solve  fi^AQ  =  u*  for  A Q.) 

7.  Perform  the  ^-inversion  for  leading  edge  cut  in  the  fine  mesh. 

This  is  region  !b'. 

8.  Perform  the  ^-inversion  for  region  1  c‘  above  the  projectile. 

9.  Perform  the  ^-inversion  for  region  fd*  below  the  projectile. 

10.  Perform  the  ^-inversion  for  region  ’ e/  the  trailing  edge  region. 

The  remaining  portion  of  this  section  serves  as  an  explanation  of  the  ten 
steps  involved  in  taking  one  integration  step  on  this  modular  grid.  Steps 
1  and  2  are  self-explanatory  except  to  say  that  line  ABG  for  Step  1  and 
line  ABF  for  Step  2  are  excluded  from  the  ^-inversion  at  this  point.  Step  3 
involves  the  ^-inversion  around  the  projectile  including  the  outer  boundary 
of  Region  1  (GBF) .  Note  that  when  the  intermesh  boundaries  are  involved, 
a  reference  side  must  be  selected  for  computing  the  geometric  data  (metrics, 
Jacobian,  etc.)  and,  in  this  case,  the  reference  side  is  taken  to  be  the 
Region  l'side.  Step  4  does  a  ^-inversion  along  the  leading  edge  line  from 
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the  inflow  boundary  to  the  projectile.  Since  in  the  fine  mesh  this  com¬ 
putational  direction  is  actually  the  negative  C  direction,  the  metrics— 

^x'  ^y '  must  used  for  (?x,  during  this  inversion  when  working 

in  Region  1.  Note  that  for  this  inversion  there  is  a  sudden  change  in  the 
grid.  As  a  result,  at  this  interface  point  between  the  fine  and  coarse 
mesh,  special  differencing  is  required  as  discussed  in  relation  to  the  model 
problem  presented  earlier.  This  is  further  discussed  in  the  next  section. 

Step  5  does  the  ^-inversion  along  the  trailing  edge  cut.  This,  of  course, 
involves  the  evaluation  of  the  RHS  along  the  intermesh  boundary  which  requires 
special  difference  operators  in  the  £ -direction.  At  this  point,  all 
^-inversions  are  complete.  Step  6  performs  a  vertical  inversion  in  region 
’a’  taking  vertically  up  to  be  the  positive  computational  coordinate.  Since 
this  is  opposite  to  the  actual  £— direction  in  Region  3  the  £— metrics  in 
Region  3  must  be  scaled  by  -1  for  this  inversion.  Note  that  again  this 
sweep  is  through  an  intermesh  boundary.  Following  Step  6,  A Q  is  known  in  • 
region  ’ a /  which  includes  point  L.  This  value  of  AQ  is  required 
for  Step  7  which  performs  the  ^-inversion  for  the  leading  edge  cut  in  the 
fine  mesh  only.  This  AQ  serves  as  a  known  boundary  condition  for  the  £-inver- 
sion  of  Step  7.  The  result  of  Step  7  is  to  generate  Aq  along  the  leading 
edge  cut  in  the  fine  mesh  (region  'b')  .  Step  8  performs  a  ^-inversion  from 
body  to  outer  boundary  in  region  *c'  above  the  projectile.  Again  an  inter¬ 
mesh  boundary  is  crossed  in  this  inversion.  Step  9  performs  a  ^-inversion 
from  body  to  projectile  centerline  (region  'd').  This  step  also  involves 
an  intermesh  boundary.  Step  10  involves  a  £— inversion  from  the  projectile 
centerline  vertically  up  to  the  outer  boundary  of  Region  2  (region  'e'). 
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Note  that  this  inversion  crosses  three  intermesh  boundaries.  The 
inversion  direction  is  opposite  of  the  actual  ^-direction  in  Region  3  and 
the  lower  part  of  Region  1  so  the  metrics  again  must  be  scaled  by  -1  in 
these  regions.  Following  these  10  steps  and  the  application  of  boundary  con¬ 
ditions  at  the  projectile  surface  and  the  outflow  boundaries,  the  value  of 
Aq  is  known  everywhere.  Thus,  an  integration  step  can  be  completed. 

The  next  section  deals  with  the  actual  intermesh  operators  employed. 


V.  INTERMESH  OPERATORS 

This  section  details  the  approach  for  performing  an  inversion  across 
an  intermesh  boundary.  Special  differencing  is  required  for  convective 
terms,  viscous  terms  and  smoothing  terms.  Each  of  these  is  discussed  in  a 
general  setting. 

A.  Convective  Terms 

A  typical  convective  term  for  the  strong-conservation-law  form  of  the 
governing  equations  can  be  expressed  as 
!G*)a 

where  a  represents  one  of  the  computational  coordinates,  G  represents  geometry 
related  information  such  as  combinations  of  metrics  and  Jacobian,  and  <J>  re¬ 
presents  the  physical  flow  quantity.  Consider  the  application  of  such  a  term 
at  an  intermesh  point,  p,  as  shown  in  Fig.  9.  Express  this  term  as 

<s+>a  -  «„♦ +  “„♦* 

where  6  is  the  arc  length  along  the  coordinate  line .  In  general  <j>  is 

continuous  across  the  interface  but  G  is  not.  As  a  result,  G  and  4  are 

a  a 

evaluated  with  one-sided  finite  differences  on  whichever  side  is  taken  to 
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be  the  reference.  This  choice  is  arbitrary  so  long  as  the  same  choice  is 

made  for  all  terms  in  the  equation.  <p  is  evaluated  at  point  p  and  <j>  is 

4 

approximated  by  the  second-order  accurate  formula 


<p  -  a  <p  +b<f>  +  c<f>  , 

4  p-1  Tp  Yp+l 


where 


(5) 


V  -  A 


c  = 


A(A  +  V) 


V(A  +  V)  '  D - AV  ' 

Thus  the  result  for  the  convective  term  is 

(G*’o  -  (GV’Vl  *  (GC.  +  “o’31  *p  +  («0=»Vl 

which  retains  tridiagonal  structure . 

B.  Viscous  Terms 

A  typical  viscous  term  may  be  generally  expressed  as 

(U  G  <t>  ) 

p  ya  a 

where  G,<f>  are  as  before  and  p  represents  a  viscosity  or  conductivity  factor. 
This  expression  can  be  transformed  to  an  arc  length  derivative  as 


<V  s  ♦„><.  "  [VaG  *  woGa  +  "“acA  * 


(6) 


The  factors  4^,  G^,  4^^  are  all  evaluated  with  one-sided  finite  differences. 

and  (p^  are  approximated  by  the  same  formula  used  for  the  convective  term 
example.  is  approximated  by  the  second-order  accurate  formula. 


~  d  K  i  +  e  $  +  f  <f>  ,, 

44  p-1  p  p+1 

where 

d  =  V (A  +  V)  '  e  =  "  w  ' 

The  result  of  these  one-sided  formulas  gives 


(7) 


f  = 


A(A  +  V) 
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(U  G  Vo-  I(&soc‘)l,p-l  +  <Giob  +  Vo  *  “oo'^p  +  <Gaoc)Vl)(*Vl  +  b*p  +  V 


+  (yG4  )  (d<J>  +  e<j>  +  f<J)  ) 

U  p— 1  p  p+1 


which  again  is  consistent  with  the  required  tridiagonal  structure.  All  of 
the  viscous  terms  have  this  form  and  are  treated  in  this  fashion  at  the 
intermesh  boundaries. 

C.  Smoothing  Terms 

The  Beam-Warming  algorithm  does  not  operate  effectively  without  added 
smoothing  terms.  The  needed  presence  of  these  terms  requires  interface 
operators  for  them,  too.  On  the  implicit  or  left-hand  side  of  the  equations, 
a  second-order  smoothing  term  is  employed.  This  term  is  identical  in  form 
to  the  typical  viscous  term,  Eq.  (6),  with  ]i  replaced  by  1.  Thus,  no  further 
discussion  of  this  term  is  needed.  On  the  right-hand  side  of  the  equations, 
a  fourth-order  smoothing  term  is  employed.  On  the  interior  of  a  region 
(i.e.,  at  point  i  =  j  or  k  in  Fig.  9),  this  term  has  the  form 


to4  r?  -  *i+2  +  *1-2  -  4*i+i  -  4*i-i +  6*i 

3  Ot 

where  <f>  is  the  conservative  dependent  variable  vector  unsealed  by  the 
transformation  Jacobian.  This  term  may  be  expressed  as 


(8) 


4  3  2 

Aa  — —  (cf>  )  =  Aa  (4> 

.  2  Yaa  aa1(1 

3a  i+I 


+  <J)  —2d)  ) 

Yaa.  .  Yaa. 

i-l  i 


where 


V,,  ■  *  *4-1  ‘  2*l 


when  t  =  i+1,  i-l,  and  i.  Although  this  is  perfectly  equivalent  to  Eq.  (8) 
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when  applied  to  points  other  than  p-1,  p,  and  p+1,  this  formula  provides 
a  mechanism  for  making  a  consistent  fourth  order  smoothing  across  the 
interface  boundary.  For  example,  if  applied  to  point  p-1,  the  result  is 


Act2  (<)) 


aa 


+  <t>, 


aa 


~  2  </> 


p-2 


aa 


p-l 


where 


=  T  +  ^  1  '2i9 

aotp_2  pj  p-i  p-2 


2 

Aa  <b  =  <b  +  <t>  -  2<b 

vaa  ^  yp-2  p  p-l 


and  (|)  is  determined  just  as  in  Eq.  (6)  with  y  =  G  =  1  with  the  result 
P 

2 

(b  =  4  (b  +  4  d) 

raa  aaY4  aY44 

P 

where  (p ^  and  <p^  are  differenced  according  to  Eqs.  (5)  and  (7)f  respectively. 

Note  that  4  and  4  must  be  evaluated  with  one-sided  differences  on  the 

aa  a 

left  of  the  interface  in  this  case.  An  identical  approach  is  used  to 
evaluate  the  smoothing  term  at  point  p+1.  That  is,  use 

Aa2  (<}>  +<})  -  2<p  ) 

aa  „  Taa  Taa  , 

p+2  p  p+1 


where 


2 

Aa  4>  =4>  „  +  4>  ,-24>  „ 

aap+2  p+3  P+1  p+2 


ll2*«  ■  Vs +  *p '  2Vi 


and  <p  is  determined  as  before  except  that  4  and  4  must  be  evaluated 

aa  aa  a 

P 

on  the  right  of  the  interface. 

The  only  remaining  task  is  to  determine  the  smoothing  term  applicable 
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at  point  p.  This  is  done  by  satisfying  a  global  conservative  property 
(see  Ref.  4)  along  the  typical  grid  line  indicated  in  Fig.  9.  This 
requires  that  the  sum  of  the  smoothing  terms  at  all  the  points  along 
this  line  equals  zero.  Let  S  represent  the  smoothing  term  of  interest. 
Then  for  t  =  p-3,.,.,p+3  becomes 


Sp-3  '  Vs  '  4V4  +  6V3  '  4*p-2  *  Vl 

Sp-2  '  *p-4  '  4V3  +  6V2  ‘  “Vl  *  *P 

Sp-1  '  *p-3  -  ^p— 2  +  Vp-1  -  SL  *p  +  Vp+1 


s  =  <p  _  +  (3  -  a  -  y  ) <p  +  (6  +  6_  -  2 ) 4>  +  (3 

p  p— 2  L  R  p“l  L  R  p 


a  -  Y  )  d)  +  6 
R  '  L  p+1  p+2 


S  =  y  <b  -fid)  +  a  < b  _  -  4<J)  ^  +  <J> 

p+1  TRyp-l  Ryp  R  p+1  Yp+2  p+3 

S  =  d>  -  4d>  ,  +  6d>  ^  -  4d)  ^  +  d> 

p+2  yp  yp+l  yp+2  yp+3  Yp+4 

S  =  cb  “  4<b  +  6 (p  —  4(p  +  (J) 

p+3  yp+l  yp+2  p+3  yp+4  Yp+5 


where 


a 


Yt 


a_ 


R 


Yt 


5  +  Aa2/^  a  +  4 2  d 
\  aaT  a 

L  i-i 


2  -  Aa  (4  b  +  4  e 
\  aa  a 

L  L 


Aa^/4  c  +  4 2  f) 
'  aaL  aL  ’ 


.2/  2 
5  +  Aa  (4  c  +  4  f 
\  aa^  ol 

x  R  B 


2  -  Aa2/ 4  b  +  42  e 

l  aa_  a_ 

'  R  R 


Aa2/ 4  a  +  42  d) 
\  aaR  aR  / 


) 

) 

) 

) 


The  constants  a,  b,  c,  d,  e,  f  are  those  used  in  Eqs .  (5)  and  (7).  The 
subscripts  L  and  R  indicate  on  which  side  of  the  interface  boundary  the 
arc  length  derivatives  are  evaluated. 
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The  smoothing  philosophy  indicated  here  is  employed  at  all  inter¬ 
face  boundary  points  and  their  immediate  neighbors. 

VI .  CONCLUDING  REMARKS 

The  tubular  projectile  code  has  been  extensively  modified  to  incor¬ 
porate  an  efficient  grid  system  consisting  of  three  distinct  modules. 

The  Beam-Warming  implicit  algorithm  has  been  modified  to  apply  to  the 
regions  of  module  interfaces.  Feasibility  of  such  module  interfaces  with 
discontinuous  metrics  was  tested  on  a  model  problem  with  positive  results. 
New  operators  were  derived  and  employed  for  use  at  interface  boundaries. 
These  operators  include  applications  for  convective  terms ,  viscous  terms , 
and  smoothing  terms.  The  topology  for  implementation  of  the  new  factored 
algorithm  is  complex  and  was  described  in  Section  IV.  The  resulting 
computer  code  is  quite  different  in  character  than  the  original  version, 
though  many  similarities  exist  in  its  structure;  its  size  is  roughly 
three  times  the  size  of  the  original.  This  code  currently  resides  on  the 
Ballistics  Research  Laboratory  Launch  and  Flight  Division's  VAX  computer. 

An  extensive  debugging  effort  still  remains,  following  which  some 
test  cases  will  be  computed  and  compared  to  existing  data. 
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Figure  1.  Hollow  Projectile 
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Figure  2.  Typical  C-Grid  for  Hollow  Projectile  Problem 
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Figure  3.  Schematic  of  New  Grid  System 
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DEPENDENT  VARIABLE, 


Figure  4.  Typical  Solution  of  Nonlinear  Viscous  Burgers'  Equation 


r 

x=o 


- <j>— 0-0-0-0-0-0-0-<j) 

X=  1 -8  X=  1 


Figure  5.  Dual  Mesh  for  Concept  Testing 
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PHYSICAL  POSITION  IN  GRID, 
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Figure  7.  Actual  Grid  System  for  a  Sharp  Leading 
and  Trailing  Edge  Hollow  Projectile 
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INFLOW  BOUNDARY 


Figure  8.  Integration  Scheme  Topology 
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Figure  9.  Arbitrary  Coordinate  Line  Through  an  Intermesh  Boundary 
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USER  EVALUATION  OF  REPORT 


Please  take  a  few  minutes  to  answer  the  questions  below;  tear  out 
this  sheet,  fold  as  indicated,  staple  or  tape  closed,  and  place 
in  the  mail.  Your  comments  will  provide  us  with  information  for 
improving  future  reports . 

1 .  BRL  Report  Number _ 

2.  Does  this  report  satisfy  a  need?  (Comment  on  purpose,  related 
project,  or  other  area  of  interest  for  which  report  will  be  used.) 


3.  How,  specifically,  is  the  report  being  used?  (Information 
source,  design  data  or  procedure,  management  procedure,  source  of 
ideas,  etc.) _ 


4.  Has  the  information  in  this  report  led  to  any  quantitative 
savings  as  far  as  man-hours/contract  dollars  saved,  operating  costs 
avoided,  efficiencies  achieved,  etc.?  If  so,  please  elaborate. 


5.  General  Comments  (Indicate  what  you  think  should  be  changed  to 
make  this  report  and  future  reports  of  this  type  more  responsive 
to  your  needs,  more  usable,  improve  readability,  etc.) _ 


6.  If  you  would  like  to  be  contacted  by  the  personnel  who  prepared 
this  report  to  raise  specific  questions  or  discuss  the  topic, 
please  fill  in  the  following  information. 

Name: 


Telephone  Number: 
Organization  Address: 


